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Abstract 

We study the rotor router model and two deterministic sandpile models. For the rotor router 
model in Z d , Levine and Peres proved that the limiting shape of the growth cluster is a sphere. 
For the other two models, only bounds in dimension 2 are known. A unified approach for these 
models with a new parameter h (the initial number of particles at each site), allows to prove a 
number of new limiting shape results in any dimension d > 1. 

For the rotor router model, the limiting shape is a sphere for all values of h. For one of the 
sandpile models, and h = 2d — 2 (the maximal value), the limiting shape is a cube. For both 
sandpile models, the limiting shape is a sphere in the limit h — ► — oo. Finally, we prove that the 
rotor router shape contains a diamond. 

1 Introduction 

In a growth model, there is a dynamical rule by which vertices of a graph are added to an initial 
collection. The existing literature on growth models deals mainly with stochastic growth models. 
Among many stochastic growth models, there are for example the Eden model, the Richardson model, 
first and last passage percolation and diffusion limited aggregation. An introduction to stochastic 
growth models and limiting shape theorems can be found in [6] . 

An example related to the models in this paper, is the internal diffusion limited aggregation 
(IDLA) model. One adds particles one by one to the origin, letting each particle perform a random 
walk that stops when it hits an empty site. The growth cluster is then the collection of sites that 
contain a particle. 

The IDLA model on Z d has been studied by Lawler, Bramson and Griffeath. In their paper from 
1992 [5], they prove that the limiting shape of the growth cluster for this model in any dimension is 
a sphere. Lawler [TO] estimated the speed of convergence. 

The three deterministic models to be discussed in this paper can be viewed as deterministic 
analogues of IDLA, and are all examples of a general height-arrow model studied in [3]. The closest 
analogue is the rotor router model, proposed by Propp (see [8]). Once more, particles are added to 
the origin of the d-dimensional lattice Z rf , but now they perform a deterministic walk as follows: at 
each site, there is an arrow present, pointing at one of the 2d neighboring sites. If a particle at this 
site finds it already occupied, then it takes a step in the direction of the arrow, and the arrow rotates 
to the next position. 
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In our other models, occupied sites hold the particles until they can be sent out in pairs, in the 
directions of a two-pointed arrow (double router model), or groups of 2d, one in every direction 
(sandpilc model). In the last model, exactly the same amount of particles is sent in every direction. 
In the other two models, the difference is at most one. 

The last two models are deterministic versions of sandpilc models. A stochastic sandpilc model, 
where addition occurs at a random site, has originally been proposed to study self-organized criticality 
PQ. Soon after, Dhar [4] introduced the abelian sandpile model, that has the symmetric toppling rule 
("toppling" of a site means that particles are sent out), which is now the most widely studied variant. 
Manna's sandpile model [16] is doubly stochastic: the addition site is random, and moreover, in a 
toppling two particles are sent to randomly chosen neighbors. 

The rotor router model on 7L d has been studied by Peres and Levine [TTJ fTS] . They have found 
that the limiting shape of this model is also a sphere, and give bounds for the rate of convergence. 
This result is at the basis of a number of new results in the present paper. Recently, a new paper of 
these authors appeared [13], further extending their results. 

The deterministic abelian sandpile model has been studied on a finite square grid [HI [TZl [20] , 
with emphasis on the dynamics, but hardly as a growth model. Le Borgne and Rossin [2] found some 
bounds for the limiting shape for d = 2. 

From the above introductory description, it should be clear that all the deterministic models we 
introduced here are closely related. In Section (2] we define the above models on Z d in a common 
framework, to allow comparisons of the models. We introduce a parameter h to parametrize the 
initial configuration. For all models, a (large) number n of particles starts at the origin, and spreads 
out in a deterministic way through topplings. The rest of the grid initially contains a number h of 
particles at each site, where h can be negative. One can imagine a negative amount of particles as 
a hole that needs to be filled up; admitting negative particle numbers will be helpful in comparing 
the different models. Once every site is stable, (i.e., has a number of particles at most some maximal 
allowed number), the growth cluster is formed. 

We present pictures, obtained by programming the models in Matlab, of growth clusters for finite 
n and d = 2. For all models, we obtain beautiful, self-similarly patterned shapes. The appeal of 
these pictures has been noted before, e.g., the sandpile pictures for h = 0,-1 can be found on the 
internet as "sandpile mandala" , but the patterns are so far unexplained. Similar patterns arc found 
in the so-called sandpile identity configuration [2l [14] . 

Finally, we explain our main limiting shape results for each model. For the sandpile model, we 
obtain that the limiting shape is a cube for h = 2d — 2, (Theorem 14. 1|) . and a sphere as h — > — oo 
(Theorem 14. 7[) . This last result also applies for the double router model. We remark that the sandpilc 
model with h — > — oo, strongly resembles the divisible sandpile introduced by Levine and Peres [13] , 
of which the limiting shape is also a sphere. 

For the rotor router model, we find that the limiting shape is a sphere for all h (Theorem 14. 3j) . 
Finally, we generalize the bounds of Le Borgne and Rossin for growth cluster of the sandpile model, to 
all h and d. As a corollary, we find that the rotor router shape contains a diamond, and is contained 
in a cube. 

The rest of the paper contains proofs of these results. In Section [3J we derive inequalities for the 
growth clusters. In Section [U we prove the various limiting shape theorems. 

2 Model definitions and results 

Before introducing each of the three models, we present some general definitions that apply to each 
model. All models are defined on Z d . We define a configuration 77 = (H, T, D) as consisting of the 
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following three components: the particle function H : Z d — > Z, the toppling function T : Z d — > Z and 
the direction function D : Z d -> P = {0, 1, . . . , 2d - 1}. We will use the 2c? possible values in P to 
indicate the 2d unit vectors eo . . . e2<j-i- The results will not depend on which value is assigned to 
which vector, as long as the assignment is fixed. 
We call a configuration allowed if 

1. For all x, T(x) > 0, 

2. For all x with T(x) > 0, if(x) > 0, 

3. For all x with T(x) = 0, if(x) > ft, 

4. For all x, D(x) e P. 

Here, /i £ Z, the "background height", is a model parameter. 

Each model starts with initial configuration rjo = r]o{n, h), which is as follows: Hq(x) = n if x = 0, 
and Hq(x) — h otherwise, Tb(x) = for all x, and -Do( x ) € P for all x. Observe that r/o is allowed. 

We now define a toppling as follows: 

Definition 2.1 A toppling of site x in configuration rj consists of the following operations: 
. T(x) -» T(x) + 1, 

• -ff (x) — > H (x) — c, with c < 2d some value specific for the model, 

• H(x + e») -> if (x + e. t ) + 1, with i = (P(x) + 1) mod 2d, ... , (-D(x) + c) mod 2d, 

• P(x) (D(x) + c) mod 2c?. 

In words, in a toppling c particles from site x move to c different neighbors of x, chosen according 
to the value of D(x) in cyclic order. We call a toppling of site x legal if after the toppling H(x) > 0. 
We call a site x stable if H (x) < h max , with h max = c — 1, so that a site can legally topple only if it 
is unstable. 

We can now define stabilization of a configuration as performing legal topplings, such that a stable 
configuration is reached, that is, a configuration where all sites are stable. We call this configuration 
the final configuration rj n = rj n {n, h), which will of course depend on the model. To ensure that the 
final configuration is reached in a finite number of topplings, we impose for each model h < h max . For 
each model, it is then known that r\ n does not depend on the order of the topplings. This is called 
the abelian property. The abelian property has been proved for centrally seeded growth models in 
general ([5], Section 4), and for the sandpile model in particular [18]. 

In the remainder of this paper, we will choose various orders of topplings. In Section [4.21 we will 
even admit illegal topplings. In Section l4.4i we will organize the topplings into discrete time steps, 
obtaining after time step t a configuration 77*, consisting of iJ*, T* and D*. 

We define two growth clusters that are formed during stabilization, denoting by x° the unit cube 
centered at x (i.e., x D = {y : y = z + x, z e [—5, |] rf }): 

Definition 2.2 The toppling cluster is the cluster of all sites that have toppled, that is, 

T n = |J x°, 

x:T„(x)>0 

and the particle cluster is the cluster of all sites that have been visited by particles from the origin, 
that is, 

v n = r n u (J x D . 

x:i? rl (x)>h 
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A cluster as a function of n has a limiting shape if. appropriately scaled, it tends to a certain shape as 
n — > oo, in some sense to be specified later. By the model definition, all clusters are path connected. 
Observe that from Definitions 12.11 and 12.21 it follows that 

T n C V„ C T n U dT n , (1) 

with dT n the exterior boundary of T n . We also define the "lattice ball": 

n 

Bn=U X ° ( 2 ) 

i=l 

where the lattice sites xi , X2 , . . . of Z d are ordered in such a way that the Euclidean distance from 
the origin is non-decreasing. 



2.1 The rotor-router (RR,/i) model 

For the rotor router model, c = 1, so that in each toppling only one particle moves to a neighbor. 
Therefore, h max = 0, and the model is defined for h < only. In words, every site holds the first 
\h\ particles that it receives, and sends every next particle to a neighbor, choosing the neighbors in 
a cyclic order. This means that after 2d topplings, every neighbor received a particle. Instead of a 
direction function with numerical values, we can think of an arrow being present at every site. In a 
toppling, the arrow is rotated to a new direction, and the particle is sent in this new direction. 

Propp, Levine and Peres studied the case h = — 1. Lcvine and Peres have proven that for h = — 1, 
the limiting shape of the rotor router model is a Euclidean sphere. More precisely, they showed f [T21, 
Thm. 1.1): 

lim A(n- 1 / d V„ A B) = 0, (3) 

n— *oo 

where A denotes d-dimcnsional Lcbesgue measure, A denotes symmetric difference, and B is the 
Euclidean sphere of unit volume centered at the origin in M. d . This result is independent of Dq, and 
of any assignment of different unit vectors to possible values of -D(x). Note that the scaling function 
is necessarily n^ 1 ^, since |V n | = n by model definition. 

In [8], a picture of D n is given for h = — 1 and n = 3 million and Dq constant. The picture shows 
a circular shape with an intriguing seemingly self-similar pattern. Curiously enough, this picture 
actually shows the shape of T n rather than V n , indicating that the limiting shape of T n is also a 
sphere. This would be a stronger statement than ([3]), but it remains as yet unproven. 

It has been noted that the shape of V n for the rotor-router model is remarkably circular, that is, 
as close to a circle as a lattice set can get, for every n. However, the shape of T n has not been studied 
before. We programmed the model for several values of h and n, and observed for all these values 
that V n \ %i is concentrated on the inner boundary of V ra . 

Our first main result of the rotor router model is the generalization of ^ to all h < — 1, stated 
in Theorem 14.31 The proof in fact uses ([3]) as main ingredient. 

By an entirely different method, in fact as a corollary of Theorem 14.81 we moreover obtain the 
result that the rotor router shape contains a diamond of radius proportional to { 2 d-i-h) 1 ^° ^ ^ s 
surprising that, in spite of the limiting shape of the RR,/i model being known, this is still a new 
result. But, in the words of Lcvine and Peres, the convergence in §5§ does not preclude the formation 
of, e.g., holes close to the origin, as long as their volume is negligible compared to n. The only 
comparable previous result is that for d = 2 and h = — 1, the particle cluster contains a disk of radius 
proportional to n 1 / 4 [IT]. 
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2.2 The abelian sandpile (SP,h) model 

For the abelian sandpile model, c = 2d. Therefore, h max = 2d — 1, and the model is defined for 
h < 2d — 1. In each toppling, one particle moves to every neighbor, so that in fact the value of -D(x) 
is irrelevant. 

It follows that for the SP model, we can specify (HJ to 

V n = T n UdT n . (4) 

The SP model as a growth model has received some attention in the cases h = — 1 ( "greedy" sandpile) 
and h = ("non-greedy"), for which pictures of H n can be found [8]. It is noted that the shape does 
not seem to be circular. In Figure [IJ we show a family of sandpile pictures for a range of values of 
h, obtained by programming the model in Matlab. We see the number of symmetry axes increasing 
as h decreases. The shape appears to become more circular as h decreases. The shape for h = 2 is 
observed to be a square, but for the other values of h it seems to tend to a more complicated shape. 

Again, we find that V n \ T n is concentrated on the inner boundary of V n , for all observed values 
of h. 

Figure 1: H n for the sandpile model with different h values. The number of particles ranges from n = 
60,000 (h = 2) to 1,000,000 (h = —20). The gray-scale colors are such that a lighter color corresponds 
to a higher value of H(x). 

For the sandpile model, we have the following results. Theorem l4. 1 1 states that indeed the toppling 
cluster for the SP,2d — 2 model is a (d-dimensional) cube, and the particle cluster tends to a cube as 
n — ► oo. 

Theorem 14.71 states that for h — > -co, the limiting shape is a sphere. 

Finally, we generalize some bounds for the scaling function that have been obtained by Le Borgne 
and Rossin [5] in the case d = 2, < h < 2, to all d and h. This result is formulated in Theorem 14.81 
Our proof moreover allows to deduce that V n for the sandpile model is simply connected for all n. 

2.3 The double router (DR,/i) model 

In the double router model, c = 2. Therefore, h max = 1, and the model is defined for h < 0. In 
a toppling of this model, two particles arc sent out in two different directions, such that after d 
topplings, every neighbor received a particle. 

Many variants of this model are possible, e.g., for d — 3 one could choose c = 3, which amounts 
to dividing the 6 neighbors in 2 groups of 3. However, to avoid confusion, we only discuss the above 
explained variant. 

Figure 2: H n for the DR model with h = (n = 110,000), h = -1 (n = 100,000), and h = -5 
(n = 400,000), and -Do( x ) = for all x. Sites with height 1 are colored white, sites with height or 
negative are colored black. 

Figure [2] shows H n for h = 0, d = 2 and n = 110.000. We ordered the unit vectors as eo = left, 
ei = right, e2 = up, e3 = down. Initially, -D(x) = for all x. Figure [3] shows D n for the same case, 
to indicate that for this model we find complex patterns both for H n and D n . 

As for the sandpile model, we see a symmetric, yet not circular shape, with notable flat edges. 
In Figure [2l also some other values of h are shown, again indicating that the shape seems to become 
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Figure 3: D n for the DR model with h = (n = 110,000), and Dq(x) = for all x. White corresponds 
to D(x) = 2, black to D(x) = 0. 

more and more circular with decreasing h. Indeed, Theorem 14. 71 states this fact for both the sandpile 
model and the double router model. 

Again, we observe that V„ \ T n is concentrated on the inner boundary of V n , for all observed 
values of h. 

3 Comparing the models 

From the model descriptions, it is clear that one SP toppling always equals d DR topplings, as well as 
2d RR topplings. Furthermore, our toppling definition (Definition 1 2. 1| ensures that one DR toppling 
always equals two consecutive RR topplings. We exploit these relations to compare the clusters for 
different models. We will with a subindex RR, DR or SP indicate the model that was used to obtain 
the growth cluster, and also use a subindex h because we compare different h values. 

Proposition 3.1 For every fixed h (for which the model is defined), d and n, 

1- Vn,SP,h C Vn^DR.h C V n ,RR,h, 

2. For all models i = RR, DR and SP: V n ,i,h-i C Vn,i,h, 

3- V n ,RR,h-l C V n ,DR,h, and V n> RRfr-{2d-l) Q Vn,SP,h- 

Proof. The proof makes use of the abelian property of all models, that is, the property that the 
stabilized configuration rj„ does not depend on the order of topplings. We are therefore free to 
choose a convenient order. Furthermore, a site x that at some instant during stabilization has either 
T(x) > 0, or T(x) = and H(x) > h, must belong to V„ in the final configuration, because by 
further topplings either H(x) or T(x), or both, increase. 

part 1. The initial configurations for all these three models are the same. We first compare V n ,sp,h 
with Vn.DR.h- We choose, for both the SP and the DR model, to first perform all legal SP-topplings. 
Since every SP-toppling consists of d DR-topplings, these are legal topplings for both models. The 
configuration is now stabilized for the SP model, but there can be sites in V n ,SP,h that are not stable 
in the DR model, since h max is 2d — 1 for the SP model and 1 for the DR model. Therefore, in the 
MD model, possibly more topplings follow, so that V n ,sp,h Q V n ,DR,h- 

We compare V n ,DR,h with Vu^rr^u in an analogous manner, this time using that h max is 1 for the 
DR model and for the RR,h model, to get V n ,DR,h Q Vn,RR,h- 

part 2. We start with comparing the RR,/i — 1 model with the RR,/i model. We choose, for both 
initial configurations, to first perform all RR-topplings that stabilize the initial configuration of the 
RR,/i — 1 model. These are legal topplings for both configurations. The configuration is now stable 
for the RR./i — 1 model, but all sites that have H (x) = for this model, have H (x) = 1 for the RR./i 
model. Therefore, for this model more topplings would follow, so that V n ,RR,h—i Q Vn,RR,h- The 
same reasoning can be applied for the DR and SP model. 

part 3. We first compare V n ,DR.h with V n rr h-i- We choose, for both models, to first perform 
all DR-topplings that would stabilize the initial configuration with height h — 1. These are legal 
topplings for both models. Then for both models, more topplings are needed. For the RR,/i — 1 
model there can be sites with H(x) = 1, that are unstable. For the DR,/i model, the same sites have 
.ff(x) = 2, and therefore are also unstable. But since one DR-toppling equals two RR-topplings, the 
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set of sites where the configuration changes by the extra topplings for the RR,/i — 1 model, is a subset 
of those for the DR,/i model. Therefore, V n ,RR,h—i f= V n ,DR,h- The argument for the sandpile model 
is similar. □ 



4 Limiting shape results 

4.1 The sandpile model in Z d , with h = 2d—2 

As noted in Section[2j Figure Q] indicates that the limiting shape for the sandpile model with h = 2 and 
d = 2 is a square. This section contains the proof of a more general statement for arbitrary dimension, 
that is, we prove that indeed the toppling cluster for the SP,2d — 2 model is a (d-dimensional) cube, 
and the particle cluster tends to a cube as n — > oo. However, we have no explicit expression for the 
scaling function f(n), so that the scaled clusters f(n)V n and f(n)T n would tend to the unit cube, 
just that this scaling satisfies n" 1 < f(n) < — |. Based on the calculations, we believe that 

f(n) is 0{n-^ d ). 

Theorem 4.1 Let <t(r) be the cube [J K . max . \ Xi \< r x °- For every n in the d-dimensional SP,2d — 2 
model, there is an r n such that 

T n = £(r„), 

and 

V n = <t(r n )Ud<t(r n ). 

For all n, this r n satisfies 

-n 1/d - - <r n <n. 
2 2 



We start by outlining the case d = 2 as an example, for the sake of clarity. When simulating the 
model for several small values of n - this can even be done by hand - one notices that the configuration 
is always such that it contains a central square with all boundary sites full, except for the corner 
sites, which have height h max — 1. Outside this square, all sites have height h max — 1. We will call 
such a rectangular boundary a critical boundary. This is all the information we need to make an 
inductive argument in n. Suppose, r\ n has a critical boundary, we add one grain to the origin and 
in the course of stabilization (to obtain r) n +i), one boundary site topples. As a consequence, all 
neighboring boundary sites topple, because they are all full and all in turn receive a grain. The two 
adjacent corner sites receive one grain, therefore become full. Therefore, after these topplings a new 
rectangular critical boundary is created. If any more boundary sites topple, we can reiterate this 
argument. We conclude that the presence of a rectangular critical boundary is stable under additions 
inside this rectangle. In the special case where additions are made only to the origin, we conclude 
by symmetry that the shape of the critical boundary will always be square. 

Below, we schematically show the argument for n = 7. 
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We start with 777 plus one extra grain at the origin. Full boundary sites are colored lightgrey, unstable 
sites grey. First the origin topples, causing the full boundary sites to become unstable. Next, these 
boundary sites topple, causing the corner sites to be unstable. The first full sites of the new boundary 
are created. 
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Next, the corner sites topple. The configuration now has a new square critical boundary. A further 
toppling of the origin, after which we obtain does not change this new boundary anymore. 

In the proof below, we generalize this critical boundary to arbitrary d and make the details more 
precise, but the idea will remain the same. 

In the course of the proof, we will need the following lemma. 



be the following configuration: H(0) = 2d, and for all other x e € r , H(x) = 
1. If the configuration <j) r is stabilized, then during stabilization, 
: max{r + 1 — max^ \xi\, 0} times. 



Lemma 4.2 Let 

hmax, otherwise H(x) = h max 
every site x topples exactly 5 Xi 

Proof. We choose to order the topplings into waves [7], that is, in each wave, we topple the origin 
once and then all other sites that become unstable, except the origin again. No site can topple more 
than once during a wave. 

In the first wave, all sites in € r topple once, because it is maximally filled. The wave stops at the 
boundary, because the sites outside £ r have at most one toppling neighbor, therefore they can not 
become unstable. All sites in £ r _i then have 2d once-toppling neighbors, so their particle number 
does not change. The sites in £ r \C r _i have less toppling neighbors, so their particle number becomes 
at most 2d — 2. 

Therefore, the effect of the first wave on <\> r , is to make it, for all sites in £ r , at most equal to 
4> r -\, with equality for all sites with £,-1. It follows that the next wave will topple all sites in £ r _i 
once. Continuing this argument, the result stated in the lemma follows. □ 



Proof of Theorem \4-l\ Wc will use induction in n. For that, wc choose to obtain the final config- 
uration as follows: the n particles are added to the origin one by one, each time first stabilizing 
the current configuration through topplings. Due to abelianness, this procedure will give the same 
final configuration as when all n particles are added simultaneously, before toppling starts. We will 
show that during this procedure, only configurations in C r are encountered, where C r is the set of 
configurations that are as follows: 

2. For all x €E T n such that max^ \xi\ = r, T„(x) = I, 

3. V n = € r Ud<Z r , 

4. For all x G V„ \ T n , H n (x) = h max =2d—l, 

5. For all x g V„, H n (x) = 2d-2. 



8 



In words, a configuration in C r is such that T n is a cube, where all inner boundary sites toppled 
exactly once, and all outer boundary sites, i.e., sites that do not belong to T n but have a neighbor 
in T n , have particle number h max . Thus, T n is surrounded by 2d maximally filled square "slabs" of 
size (2r + All sites outside V n have particle number 2d — 2, by model definition. 

After the first particle is added, H\(0) = h max , and no site has toppled yet. After the second par- 
ticle is added, the origin topples once, and all neighbors of the origin receive one particle. Therefore, 

m e Co- 

We will now show that if we assume r\ n G C r for some r, then either rj n+ \ G C r or r} n +\ G C r +\. 
To prove this, we will add a particle at the origin to configuration r\ n and start stabilizing. First we 
show that if during stabilization the configuration leaves C r , then it enters C r +i. Then we show that 
in this last case, it does not leave C r+ \ during further stabilization. 

If no site x G V n such that max^ \xi\ = r topples during stabilization, then rj n+ \ remains in C r . 
But if one such site topples, then one site in V n \ T n becomes unstable, and also topples. This site 
is in one of the 2d maximally filled slabs. If one site of such a slab topples, then the entire slab 
must topple, because all of its sites will in turn receive a particle from a toppling neighbor. After 
the entire slab toppled once, all neighbors of the slab received a particle, and the configuration can 
be described as follows: The toppled cluster is a rectangle of size (2r + l) d_1 (2r + 2), surrounded by 
maximally filled slabs, two of them cubic and the rest rectangular. 

But due to symmetry, if one slab topples once then this must happen for all slabs. We can 
choose in what order to topple the slabs; suppose we divide them in d opposite pairs. After we 
toppled the first pair of slabs, the configuration is a central rectangle of size (2r + l) d_1 (2r + 3), 
surrounded by maximally filled slabs, two of them cubic and the rest rectangular. After we toppled 
the k th pair of (meanwhile possibly rectangular) slabs, the configuration is a central rectangle of size 
(2r + l) d ~ fc (2r + 3) fc , surrounded by maximally filled slabs, so that after all slabs toppled, the toppled 
cluster is again a cube, now of size (2r + 3) rf , centered at the origin and surrounded by 2d maximally 
filled square slabs. In other words, after these topplings the configuration is in C r +i. It follows that 
if rj n G C r for some r, then rj n+ \ G C r >, for some r' > r. It now remains to show that r' can only have 
the values r or r + 1 . 

We use lemma l4~2l A configuration r\ G C r , plus an addition at the origin, has at every site at 
most the number of particles given by 4> r +i- Therefore, if we suppose r/ n G C r , then upon addition of 
a particle at the origin, by abclianness, the number of topplings for every x will at most be <5 x .r+i- 
In particular, the outer boundary sites of T n will topple at most once. 

The induction proof is now completed, therefore we now know that, for all n, r\ n G C r for some r. 
We have also shown that r„ < n. From the description of C ri we see that T n is always a cube, and 
that V n is more and more like a cube as r n increases. 

To conclude that V n tends to a cube as n — ► oo, we finally need to show that r„ increases with n. 
For this, note that in rj n every site needs to be stable. Since h = h ma x — ^, every site can accommodate 
at most one extra particle. Therefore, if r\ n G C r , then r > ^n 1 ^ — |. □ 

4.2 The rotor-router model with h < — 1 

In this section, wc will prove the following result: 

Theorem 4.3 The limiting shape of the particle cluster for the rotor-router model is a sphere for 
every h < — 1. More precisely, 
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The proof will also reveal the following about the rotor-router toppling cluster: 
Corollary 4.4 The toppling cluster of the RR,/i model contains a cluster W n , with 

^ A ( ( NTT rl/dw « AB )=°- 

Theorem 14.31 has been proven for the case h = — 1, by Levine and Peres, see ([3]). We will use 
their result to prove the theorem for other values of h. The 1 strategy of the- proof will be as follows. 
First, we introduce a slightly different version of the RR,/i model, which we will call the fc-color rotor 
router, or RR',fc model. We show that the limiting shape of this model is a sphere for all fc, with 
scaling function (r)~ • Informally, the fc-color model can be viewed as fc iterations of the RR,-1 
model. The n particles at the origin are equally divided in fc different-colored groups, and for each 
color the particles spread out until there is at most one particle of each color at each site. From ([3]), 
we then know that each color region is approximately spherical, so we get a final configuration that 
looks like fc almost completely overlapping, approximately spherical, different-colored regions. 

Then we show that V n- RR^h and V n ,RR',k, with h = — fc, differ only in a number of sites that is 
o(n). This number is at most the number of sites in V n ,RR',k where not every color is present. In the 
proof of this second point, we will stabilize rjo for the RR,/i model by first performing all topplings 
needed to stabilize the RR',— h model. Some of these topplings may be illegal for the RR,/i model, 
so we will reverse them by performing untopplings. The main difficulty in the proof is to show that 
we do reach r\ n by this procedure. First we define untopplings. 

Definition 4.5 An untoppling of site x in configuration r\ consists of the following operations: 
. T(x) -» T(x) - 1, 

• H (x) — > H(x) + c, with c according to the model, 

• D(x) (.D(x) - c) mod 2d, 

• H{x + e;) -> H{x + e;) - 1, with i = D(x), (D(x) + 1) mod 2d, ... , (D(x) + c - 1) mod 2d. 

We call an untoppling of site x legal if T(x) > before the untoppling. We see that a legal 
untoppling of site x is precisely the undoing of a toppling of site x, in the sense that the particles 
that were sent out of site x in the last toppling, return to this site, and the value of -D(x) returns to 
its previous value. 

The following proposition provides a more elaborate description of rj n in the case where we allow 
illegal topplings, that we will need in the proof of Thcorcm l4.3l fwe remark that the abelian property 
still holds for combinations of legal and illegal topplings) . We introduce the notion of optimality to 
distinguish rj n from configurations that are stable and allowed, but that cannot be obtained from 770 
by legal topplings. An example for the rotor router model is as follows. Suppose that, after reaching 
r) n by legal topplings, there is somewhere a closed loop of sites with height such that "the arrows 
form a cycle", i.e., if we topple all of them, they would all send and receive one particle. After these 
topplings (of which at least the first would be illegal), the height function is still H„, but the toppling 
number of these sites has increased by 1, so we obtained a different, stable and allowed configuration. 
We remark that a set of sites forming such a "cycle of arrows" cannot be found in T n . Since all sites 
that toppled are full, in the last toppling in this set a particle must have left the set, so that there is 
an arrow pointing out of the set (see also [T5], section III). Therefore, if there is a set of sites forming 
a cycle of arrows, then at least one of these sites did not topple. 
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Proposition 4.6 Call a stable configuration optimal if every sequence of legal untopplings leads to 
an unstable configuration. For every model, r] n is the unique optimal, stable and allowed configuration 
that can be reached from r/o by topplings, cither legal or illegal, and legal untopplings. 

Proof. We defined rj n before as the unique stable configuration that can be reached from r/o by legal 
topplings. It follows from this definition that r\ n is allowed, since if sites can only topple when they 
are unstable, then we automatically obtain for all x where T n (x) > 0, that H n (x) > 0. 

To prove that r\ n is optimal, we proceed by contradiction. Suppose that, starting from 77™, there 
is a sequence of legal untopplings such that a stable configuration £ is obtained. Then £ can be 
obtained from 770 by a sequence of legal or illegal topplings. Call T" the toppling function for this 
sequence, then < T'(x) < T„(x) for all x (because we undid some topplings of T n ), but T' ^ T n . By 
abelianness, £ depends only on 779 and T' . Thus, we can choose the order of the topplings according 
to T", to obtain £. There cannot be an order such that all topplings are legal, otherwise £ cannot be 
different from rj n . We will therefore choose the order such that first all possible legal topplings are 
performed, and then the rest (as an example, suppose T'(0) = 0. Then we must start with an illegal 
toppling, since in 770 all sites but the origin are stable). After all possible legal topplings according 
to T", we have a configuration with at least one unstable site, since we did not yet reach r\ n . In the 
remainder, this site is not toppled, because otherwise another legal toppling could have been added 
to the first legal topplings. Therefore, £ cannot be stable. 

We so far established that r/ n is optimal, stable and allowed. Now we prove that it is unique by 
deriving another contradiction. Suppose there is an other optimal, stable and allowed configuration 
£ that can be reached from 770 by topplings and legal untopplings. Call T" the toppling function for 
this sequence. Then T"(x) > T„(x) for all x. We can see this as follows: we choose the toppling 
order such that we first perform all topplings that are also in T n . If there were some sites y where 
T"(y) < T n (y), then at this point at least one of these sites would be unstable. But in the remainder 
of T", this site would not topple again, so in that case £ could not be stable. 

Call t(x) = T"'(x) — T„(x). Since ( ^ i)„, T" ^ T„. Suppose we first perform the topplings 
according to T„, and then according to r. Then we have that £ can be reached from r] n by performing 
r(x) topplings for every site x, and vice versa, that r\ n can be reached from £ by performing r(x) 
legal untopplings for every site x. But then £ cannot be optimal. □ 

Proof of Theorem \4 -3\ The RR',fc model is defined as follows: As in the KR,h model, c = 1. Initially, 
there are n particles at the origin; at every other site there are — k particles. We choose n a multiple 
of k, and divide the n particles into k groups of n/k particles, each group with a different color. The 
height H (x) of a site x is now defined as — k plus the total number of particles present at x, of either 
color. However, now we call a site stable if it contains at most one particle of each color. If a particle 
arrives at a site where its color is already present, then that particle will be sent to a neighbor in 
the subsequent toppling. Note that in this model, a site x can be unstable even if H(x) < 0. We 
furthermore restrict the order of legal topplings for this model. We say that in this model, first all 
legal topplings for the first color should be performed, then for the next, etc. Then for each color, 
the model behaves just as the RR, —1 model, with for each color a new initial T and D. Therefore, 
with this toppling order the final configuration consisting of H n , T n and D n , is well-defined. 

We now show that the limiting shape of the RR',fc model is a sphere. We perform first all the 
topplings with particles of the first color, say, red. From [12], Thm. 2.1 it follows that these particles 
will form a cluster V^/k c l° sc to the lattice ball B n /k- In fact, the number of sites in V^ k A B n / k is 
o(n). We will denote this as 

\V^ /k AB n/k \<f(n/k) = o(n). 
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D is now different from Dq. However, if next we stabilize for the blue particles, these will again form 
a cluster close to the lattice ball B n / k , since the result © does not depend on the initial D. Of 
course, this cluster V^, k need not be the same as V* , fc . 

When all sites are stable, we have V n .RR\k = Uj=i ^n/k- We a ^ so define a cluster W n ,k = 
(li=i Kx/fc- F rom the above, the number of sites in the difference of both these clusters with is 
at most kf(n/k). Thus 

\V n , RR ,, k A B n/k \ < kf(n/k) = o(n). (5) 

We will call X n = V n .RR' ,k \ W n ,fc, so that X n contains \X n \ < 2kf(n/k) sites. Sites in W n ,k contain 
k particles, sites in X n contain less than k particles. 

Now we compare this model with the RR,/i model, with h = —k. Disregarding the colors, the 
initial configuration is the same for both models. Suppose we perform in the RR,/i model all the 
topplings that are performed as above in the RR',— h model. The configuration is then as follows: 
for all x e W„, k , H(x) = and T(x) > 0, and for all x' e X n , -k < H(x') < and T(x') > 0. 

This configuration is possibly not rj n for the RR,/i model, since it is possibly not allowed, if 
there are sites x' in X n with T(x') > 0. We can reach an allowed configuration by performing legal 
untopplings; it will appear that we then in fact reach rj n . 

Suppose first that only one untoppling is needed. It might be that the neighbor y that returned 
a particle, now has less than k particles, so that H(y) < 0, while T(y) > 0. In that case, y would 
now also have to untopple. We call this process an untoppling avalanche. The avalanche stops if an 
allowed (stable) configuration is encountered, that is, if the last particle came from a site that did not 
topple. An untoppling avalanche consists of untoppling neighbors, each one passing a particle to the 
previous one. Therefore, an untoppling avalanche of arbitrary length, changes the particle number 
of only two sites in the configuration. 

In our case, at most (k — l)|<%Vt| untoppling avalanches are required. Each avalanche stops if 
the last particle came from a site that did not topple. This will change the configuration in at most 
2{k — l)\X n \ sites. After all these untopplings, the configuration is allowed and stable, i.e., we have for 
all sites with T(x) > 0, that H(x) = 0. We can only perform legal untopplings at sites with T(x) > 0, 
so that we cannot perform a sequence of legal untopplings in a closed loop of sites. Therefore, any 
further sequence of legal untopplings would increase H(x) for at least one of these sites, and thus 
lead to an unstable configuration. 

We have thus reached an optimal configuration, so we have reached r\ n . Therefore, at this point 
we can conclude 

IVn.flH.h A V n , RW<k \ < 2(k - l)\X n \ < 4k(k - l)f(n/k) - o{n). 
Combined with ([5]), this leads to Theorem 14.31 □ 

To prove the corollary, we define W n h — Ux xev rr h h (x)=\h\ x " above argument, we have, 

as for W n ,k, that \W n ,h A B n /\ h \ \ = o(n). 

Suppose, in the RR,ft. model, we first perform all topplings that are needed to stabilize for the 
KR,h — 1 model, as in the proof of Proposition 13 . 1 1 part 2. The particle cluster then contains a cluster 
Wn.h-1, with li m?woo A ((p^-V^w,, Afl)=0. 

Because the sites in W n ,h-i contain \h — 1| = \h\ + 1 particles, every site in W TCj /i_i is unstable in 
the RR,/i model. Therefore, W n ,h-l £ %i,RR.h- d 

4.3 The double router model, and the sandpile model with h — > — oo 

In discussing Figures [T] and we observed that the DR and SP shapes seem to become more circular 
as h decreases. Proposition 13 . II and Theorem 14 . 31 can indeed be combined to give the following result: 
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Theorem 4.7 The limiting shape of the SP,ft and the DR,/i model, for ft — > — oo, is a sphere. More 
precisely, 

lim limsupA ( (^r)~ 1/d V„ A B ) = 0. 



Proof. From Proposition ^. 1) for all ft. < — 1, 

V„,iifl l / l _(2d-l) C V n ,SP,h £ Vn.Rfl^, (6) 

and 

Vn^ij^-i C V n ,DR,h C Vn.jfcR,/,. (7) 

We will discuss the DR,ft, model first. 

From Theorem 14.31 we know that the particle cluster of the RR,ft model, scaled by (m) -1 i 
tends for every fixed value of ft to the unit volume sphere as n — » oo, i.e., 



lim A ( n r)- 1/d V„,^. h AB = 0, 
and thus we also know (note that since h is negative, \h — 1| = \h\ + 1) 



Jim A^J-^Kw-iAB 
Because of (O, it follows that 



1 ^ +i J 



AS = 1 



\h\ 



\h\ + l 



< limsupA ( {^)- l l d V n , DR , h Ab) <1 



\h\> "n,un,n~~j^- ^ + j , 



so that 



lim limsupA fer^K,^ AS ) = 0. 

h— >-oo n— >oo \ |«| / 



The argument is similar for the SP,ft model. □ 
4.4 The toppling cluster for the sandpile model 

In this subsection, we generalize Theorem 4 of Le Borgne and Rossin [2] , who studied the particle 
cluster of the SP,ft model for d = 2 and ft = 0, 1 and 2, to arbitrary d and ft.. 

Theorem 4.8 Let 3)(r) be the diamond Ux:£*_, |i ! |<r xD ' and C ( r ) tnc cubc Ux:max» l^o* - 

1. In the SP,ft model, for every n and ft < 2<i — 2, there is an r„ such that 

D(r„ -1)CT„C C( r „). 

and 

S)(r„) CV„c C( r „ + 1). 

2. This r„ satisfies 

i/d 

— 3 < 2r n for all n and ft < 2d — 2, 



2d-l-7i 



2r„ < ( 3=5 I + o{n l / d ) for n -> oo and ft < rf. 
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Corollary 4.9 V n ,RR,h contains the diamond D ( | y ->d-\-h j — I ) f° r au h < 0, and is contained 
in the cube £(^(^) 1/d ) for all h < -d. 

The corollary follows from combining the theorem with Proposition 13. II part 1 and 3. 

Remark 4.10 Note that the first inequality for r n agrees with that in Theorem 14. 11 for the sandpile 
model with h = 2d — 2. We can not use the second inequality in this case, because that is only valid 
for h < d. From combining Theorem 14. II with Proposition 13. II part 3, we find that for — d < h < 0, 
Vn,RR,h is contained in the cube £(n +1). 

Proof of Theorem \4-8[ part (1) We prove the inequality for T n . The inequality for V n then follows 
from (g]). 

We follow the method of Le Borgne and Rossin, who introduce the following stabilization pro- 
cedure: Starting at t = in the initial configuration with n particles at the origin, each time step 
every unstable site of the current configuration is stabilized. For example, to obtain rj 1 we topple 
only the origin, as many times as legally possible. At t = 2, we topple only the neighbors of the 
origin, et cetera. We will call this the Le Borgnc-Rossin procedure. Thus, rf contains stable sites 
that toppled at time t, and unstable sites that received grains. T* contains for every site the total 
number of topplings up to time t. For every site x, we then have 

T t+1 (x) = max{ 

where y ~ x means that y is a neighbor of x. We remark that this formula should correspond to 
formula (2) of [2], but, apparently due to a misprint, there the term H°(x) was omitted. 

If we divide the grid Z d in two ( "checkered" ) subgrids {x : J2i=i x i ^ s even} and {x : J2i=i x i 1S °dd}, 
then at every time t, at most one subgrid contains unstable sites, because the neighbors of sites of 
one subgrid arc all in the other subgrid, and we started with only one unstable site. We will consider 
next-nearest neighbor pairs x and z, i.e., x and z are in the same subgrid, at most 2 coordinates 
differ, and \\xi\ — \zi\\ = 2, or in other words, z is a neighbor of a neighbor of x, x itself being 
excluded. 

First, we will prove by induction in t that for every next-nearest neighbor pair x and z, with 
<i(x) < d(z), where d(x) is the Euclidean distance of x to the origin, T'(x) > T'(z). We remark that 
for a next-nearest neighbor pair x' and z' with d(x') = d(z'), we must have that x' is equal to z' up 
to a permutation of coordinates. The model should be invariant under such a permutation, therefore, 
in that case we have T*(x') = T'(z'). 

The statement is true at t = 0, since at t = 0, we have T°(x) = for all x. It is also true at t = 1, 
since at t = 1, only the origin topples. 

Now suppose that the statement is true at time t, and let Xi and X2 be two next-nearest neighbor 
sites that do not topple at time t. Suppose d(xi) < d(x2). Then by assumption, X) yi ~ X i ^*(yi) — 
Ey,~x, ^'(yz), because the neighbors of xi and X2 can be grouped into next-nearest neighbor 
pairs with d{y\) < d{y2)- Furthermore, we always have if°(xi) > £f°(x 2 ), since H°(0) = n. 
and H°(x) = hifx^0. 

Inserting this in ©, we obtain T t+1 (x!) > T t+1 (x 2 ), so that the statement remains true at t + 1 
for one subgrid. The other subgrid must contain only sites that do not topple at t + 1, so that for all 
x in this subgrid, T t+1 (x) = T*(x). Therefore, the statement remains true at t+ 1 for both subgrids. 

Now let x be a site with maximal r = max^ \xi\ where T(x) > (by symmetry, x cannot be 
unique). Then no sites outside £(r) can have toppled. Furthermore, all next-nearest neighbors of x 



id K« 



,0}, 



(8) 
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that are closer to the origin have also toppled, and subsequently all next-nearest neighbors of those 
next-nearest neighbors, that are still closer to the origin, etc. Then we use symmetry to find that all 
sites in 2)(r) and in the same subgrid as x must also have toppled. An example for d = 2 is depicted 
in Figure l4~4l 
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Figure \4-4l Let x be a site with maximal r = max^ |a;j| that toppled. In this example, r = 3. By 
the fact that all next-nearest neighbors that are close to the origin have also toppled, we find that sites 
y, and the origin, have also toppled. By symmetry, we find that sites z have also toppled. Then we 
conclude that all sites in S)(3) (colored lightgrey) and in the same subgrid as x have toppled, and no 
sites outside £(3) in the same subgrid can have toppled. 

Furthermore, for at least one of the neighbors y of x, we must have T(y) > 0, because T n is path 
connected. This neighbor has max^ \yt \ > r — 1. For this neighbor, we can make the same observation. 
It follows that all sites in 35 (r — 1) have toppled, so that the first part of the theorem follows. □ 

To prove the inequalities for r n , we need the following lemma: 

Lemma 4.11 Let p n be the average value of H n (x) in T n , o~ n the number of sites and (3 n the number 
of internal bonds in T n . Then 

— < Pn < 2d -I. 
On 

Proof. It suffices to prove that the configuration restricted to T n is recurrent [4j [18] . This can be 
checked in a nonambiguous way with the burning algorithm [3] . The inequality then follows from the 
fact that in a recurrent, stable configuration restricted to some set, the total number of particles in 
the set is at least the number of internal bonds in the set, and at most h max times the number of sites 
in the set. In the case of an unstable configuration, the term "recurrent" is somewhat inappropriate, 
but we still define it to indicate configurations that pass the burning algorithm. 
We use the following properties of recurrent configurations: 

1. recurrence of a configuration restricted to a set is conserved under addition of a particle to the 
set, and under a legal toppling in the set. 

2. if the configuration restricted to a set is recurrent, and we add an unstable site x to the set, 
then the configuration is recurrent restricted to the extended set. 

We will use induction in n. We will show that if rj n restricted to T n is recurrent, then r) n +i 
restricted to 7~ n+ i is recurrent. For a starting point of the induction, we choose n' = 2d — h, so that 
\T n ' \ = 1, because a configuration restricted to a single site is always recurrent. 

We can obtain by starting from r] n , adding a particle to the origin, and stabilizing. If we 
first stabilize restricted to T n , then by the first property of recurrence, the configuration restricted to 
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T n is still recurrent. If T n = T n +i, then we are now done. But it is also possible that at this point, 
there are unstable sites outside T n . 

By both properties of recurrence, if we add one of these sites to T n and stabilize with respect to 
the new set, the configuration restricted to the new set is still recurrent. Wc can repeat this step until 
there are no more unstable sites to be added. Then the new set is T n +\, and the new configuration 
is r) n+1 . □ 

Proof of Theorem \4-8\ part (2) To calculate the inequalities for r n , we observe that the minimum 
r n would be obtained if V„ would equal the cube £(r n + 1), with each site containing the maximal 
number of particles. The number of sites in this cube is (2r n + 3) d , and the number of particles per 
site is then 2c? — 1 — h, so that 

n< (2r„ + 3) d (2d - 1 - h), 

for every n. 

The maximal r n would be obtained if V„ would equal the diamond 2)(r„), containing the minimum 
number of particles. As we deduced that the configuration restricted to T n is recurrent, we calculate 
the minimum average value of H(x) on D(r n — 1), using Lemma 14.111 The number of internal bonds 
in D(r n - 1) is (2r„ - 2) d . The number of sites in £)(r„ - 1) is ^(2r„ - l) d + o(r d ) as r„ -> oo. Thus, 

n > (2r„ - 2) d - -(2r« - l) d + o{r d ), r n -> oo. 

As in the limit n->oowe also have r n — > oo, this leads to the second inequality for r n in the theorem. 

□ 

The proof of Theorem 14.81 allows to prove an interesting characteristic of V„: 

Proposition 4.12 For the sandpile model, V n is simply connected for all n. 

Proof. We first prove that T n is simply connected, by contradiction. 

Suppose that Tn.SP.h is not simply connected, that is, T n ,sp.h contains holes. This means there 
must be a site y in T n ,SP,h 'beyond a hole', that is, a site y with T(y) > 0, and at least one neighbor 
x of y with T(x) = 0, that is closer to the origin than y itself and all neighbors x' of y for which 
T(x') > 0. At least one neighbor x' of y must have T(x') > 0, because T n is path connected. Thus, 
among the neighbors of y, there must be a next-nearest neighbor pair z and z', with d(z') > d(z), 
T(z) = and T(z') > 0. But this contradicts the above derived property that for all t, for every 
next-nearest neighbor z and z' with d(z) > d(z'), T*(z) < T*(z'), as this property should also hold 
for the final configuration. 

The above does not suffice to conclude that V„ is also simply connected. Since V n = T n U 5T n , 
we must show that T n does not contain so-called fjords, i.e., places where the boundary of T n nearly 
touches itself, such that T n U 5T n would contain a hole. But when one supposes that T n does contain 
a fjord, then one can derive the same contradiction as above. Therefore, V n is simply connected. □ 

References 

[1] P. Bak, K. Tang and K. Wiesenfeld (1988) Self-organized criticality, Phys. Rev. A 38, 364-374. 

[2] Y. le Borgne and D. Rossin (2002) On the identity of the sandpile group, Discrete Math. 256, 
775-790. 

[3] A. Dartois and D. Rossin (2004) Height Arrow Model, Formal Power Series and Algebraic Com- 
binatorics, 16th conference, Vancouver. 



16 



[4] D. Dhar (1990) Self-Organized Critical State of Sandpile Automaton Models, Phys. Rev. Lett. 
64(14), 1613-1616. 

[5] P. Diaconis and W. Fulton (1991) A growth model a game, an algebra, Lagrange inversion, and 
characteristic classes, Rend. Sem. Mat. Univ. Pol. Torino 49, 95-119. 

[6] R. Durrctt (1988) Lecture notes on particle systems and percolation. Wadsworth and Brooks/Cole, 
PacAfic Grove, Calif. 

[7] E.V. Ivashkevich and V.B. Priezzhev(1998) Introduction to the sandpile model, Physica A 254, 
97-116. 

[8] M. Kleber (2005) Goldbug variations, Math. Intelligencer 27 (1), 55-63. 

[9] G.F. Lawler, M. Bramson and D. Griffcath (1992) Internal diffusion limited aggregation, Ann. 
Probab. 20(4), 2117-2140. 

[10] G.F. Lawler (1995) Subdiffusivc fluctuations for internal diffusion limited aggregation, Ann. 
Probab. 23(1), 71-86. 

[11] L. Lcvinc (2002) The rotor-router model, Harvard University senior thesis. 

[12] L. Levine, Y. Peres (2005) Spherical Asymptotics for the Rotor-Router Model in Z d , to appear 
in Indiana University Math Journal. 

[13] L. Levine, Y. Peres (2007) Strong Spherical Asymptotics for Rotor-Router Aggregation and the 
Divisible Sandpile, Preprint available at |http://www.arxiv.org/abs/0704 .0688, 

[14] S.H. Liu, T. Kaplan and L.J. Gray (1990) Geometry and dynamics of deterministic sand piles, 
Phys. Rev. A 42(6), 3207-3212. 

[15] S. Liibeck, N. Rajewsky and D.E. Wolf (2000) A deterministic sandpile revisited, Eur. Phys. 
Journ. B 13, 715-721. 

[16] S.S. Manna (1991) Two-state model of self-organized criticality, J. Phys. A: Math. Gen. 24, 
L363-369. 

[17] M. Markosovaand P. Markos (1992) Analytic calculation of the attractor periods of deterministic 
sandpilcs, Phys. Rev. A 46(6), 3531-3534. 

[18] R. Meester, F. Redig and D. Znamcnski (2001) The abelian sandpile model, a mathematical 
introduction, Markov Proc. Rel. Fields 7, 509-523. 

[19] A.M. Povolotsky, V.B. Priczzhcv and R.R. Shchcrbakov (1998) Dynamics of Eulerian walkers, 
Phys Rev. E 58, 5449-5454. 

[20] K. Wiesenfeld, J. Thciler and B. McNamara(1990) Self-organized Criticality in a Deterministic 
Automaton, Phys. Rev. A 65(8), 949-952. 



17 



This figure "mannaD.gif" is available in "gif" format from: 



http://arXiv.org/ps/math/0702450v2 



This figure "mannashapes.gif" is available in "gif" format from: 



http://arXiv.org/ps/math/0702450v2 



This figure "sandpileshapesanders.gif" is available in "gif" format from: 



http://arXiv.org/ps/math/0702450v2 



